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Abstract 

We analyze the collective surface diffusion coefficient, D c , near first-order phase transitions 
at which two low temperature phases coexist, using a lattice-gas model on a triangular lattice. 
Assuming the local equilibrium approximation, we obtain, both numerically and analytically, the 
dependence of D c on the coverage and system size near such transitions. Four specific examples 
of transitions are considered for illustration. It is pointed out that the derived results are in fact 
rather general and not restricted only to the studied model. 
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I. INTRODUCTION 



The collective (or chemical) surface diffusion coefficient, D c , is defined via the Fick's 
first law and represents a relevant transport coefficient for surface diffusion. Theoretical 
studies of D c and of the influence of lateral interparticle interactions on D c have employed 
various approaches, such as the mean-field, real-space renormalization group, and computer 
simulation techniques 

n a. 

One of the intriguing problems that has attracted particular 
attention is the effect of a phase transition on surface diffusion. Since lattice-gas models 
can be used to model the transition, they have provided a convenient framework in this 
regard [3- 14] . However, below critical temperatures ordered phases may arise due to lateral 
interactions, and sophisticated arguments should be applied to analyze surface diffusion jfj]. 

In this paper we investigate the coverage dependence of the diffusion coefficient D c at low 
temperatures for a model on a regular triangular lattice in which phase transitions between 
ordered phases occur. To this end, we consider an approximate expression for D c valid in 
the local equilibrium limit js, 7], 

and use that the quantities in this expression can be obtained from the model finite-size free 
energy, /, that we derived in a previous paper 15[ (here D® is the diffusion coefficient of 
non-interacting particles, /3 = 1/ksT is the inverse temperature, \x is the chemical potential, 
and x is the isothermal susceptibility). In the approximation of Eq. (JT]) particles are assumed 
to migrate on the surface by jumps only to nearest-neighbor sites without any correlation 
between the jumps. 

Several forms have been introduced for the configuration factor P. Nevertheless, many of 
them have an identical structure: they can be written as a sum of derivatives of / p, Q, 17]. 
Therefore, we shall consider a very simple case in which P is identified with the probability 
that a lattice bond is vacant; then, on the triangular lattice, one has 

P = l-20 + ~£, (2) 

where 9 is the surface coverage and £ is the statistical average number of occupied bonds 
per site. We will eventually show that our results do not pertain only to form (]2]) but are 
in fact general and can be readily extended to other forms of P. 
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The paper has the following structure. The studied model is introduced and its low 
temperature phases and free energy / are described in Sec. [Ill The coverage dependence of D c 
is analyzed in Sec. 1111} focusing on the derivation of finite-size formulas for the dependence. 
Extensions of the obtained results to the more general forms of P are described in a final 
section. 



II. THE MODEL 



We assume that particles can be adsorbed on a solid surface only at sites forming a regular 
triangular lattice. The system contains a rectangular array with a large but finite number, N, 
of adsorption sites. Periodic boundary conditions are applied so that the array forms a finite 
torus whose cell is given, for example, by the vectors (3n, 0) and (0, 2y / 3^) with n = 1, 2, . . . 
[the elementary lattice vectors are (1,0) and (1/2, a/3/2)]; thus, N = (3n) x (An). Each 
lattice site is either vacant or occupied by a particle. The model Hamiltonian is 6|, |7| 

H = e b N b + e t N t - fxN s , (3) 

where N b , N t , and N s is the number of occupied bonds, elementary triangles, and sites, 
respectively, and e b and Et are bond and triangle interaction energies. This model was 
already used to study surface diffusion on a triangular lattice at high temperatures for the 
special cases when e b = (with (3\e t \ below 4.8) and e t = (with (3\e b \ below 10) 
Here we consider the general case when both the bond and triangle interactions e b and e t 
are present and temperatures are sufficiently low. 

As we proved in Ref. Q, model ([3]) has four ground states [see Fig. [1(a)]: a fully vacant 
state, do, a fully occupied state, a±, and two threefold degenerate states, <ti/ 3 and 02/3, with 
the partial coverages of 1/3 and 2/3, respectively. The ground-state diagram is shown in 
Fig. [Jjb)-(d) and can be easily constructed by comparing the ground-state energies e = 0, 
e± = 3e b + 2e t — /i, ei/3 = — /x/3, and e2/3 = e b — 2/x/3. On the lines separating the ground- 
state regions of <To and <Ti, <to and 02/3, an d c\ and 0i/ 3 only the two associated ground states 
coexist, while on the remaining lines (as well as at the points where three or all four ground- 
state regions meet) there is an infinite number of ground states, yielding in fact a residual 
entropy. Each ground state a a , a = 0,1/3,2/3,1, gives rise to a unique low-temperature 
phase, p a , whose typical configuration looks clS Si SGcl ' of the ground state a a with isolated 
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FIG. 1: (a) The ground states of the model. Circles (disks) represent (occupied) sites, (b)-(d) 
The ground-state diagram in dependence on the sign (attractivity or repulsivity) of the triangle 
interaction Et- Solid lines and disks correspond to those parts where infinitely many ground states 
exist, while dashed lines to the parts with a finite number of ground states. 



'islands' of non-ground-state configurations, thus resembling the structure of a a . So, phase 
Po (pi) is fully vacant (fully occupied), while phase py 3 (P2/3) has the occupancy of 1/3 (2/3). 
This holds only if the number of ground states remains finite, i.e., within each ground-state 
region and on the lines where only two ground states coexist. Consequently, a first-order 
phase transition can take place only between phases po and pi, po and P2/3, and pi/3 and p\. 

Let us consider one of these pairs of phases to be denoted as p- and p + (specific pairs of 
phases will be considered in Sec. Ill II when applying the results to examples); the associated 
ground states are a_ and cr + , respectively. At low temperatures and near the transition point 
lit (for \fi— fi t \ < const / /3v~N), the finite-size partition function, Z, can be expressed as a sum 
of two single-phase partition functions. Each of them contains only the contributions from 
the ground state <r_ or a + and its thermal perturbations (the latter can be written in the form 
of a cluster expansion). As a result, the finite-size specific free energy, / = — (1//3N) In Z, 
reads Q, Q 

/(^ = __L ]n[ v _e-M-V)N + v +e -M+V)N] + r , (4) 



0N 
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where the degeneracy v± of phase p± is Vq = V\ = 1 and 1/1/3 — ^2/3 — 3 and the error 
r = 0[exp(— const /3\/N)} [the symbol O(x) is used to denote a term that can be bounded 
by const a;]. Moreover, the single-phase specific free energies f± are essentially equal to 
the specific energies, e±, of ground state a± because the contributions from the thermal 
perturbations of cr± are exponentially suppressed in /3 (the Peierls condition). Taking into 
account only one-site perturbations (which represent the leading corrections), one has l^] 

/a « e a - i ln[(l + e^ A ^)«(i + e -/»AflS ) i-«] f (5) 

where 

A# « = -AflJ /s = -/i, 
AH; /3 = -AHZ /3 = 3e b -fji, (6) 
A.£/"'/3 = -AF° = 6e fe + 6e t - // 
are the energy excesses of one-site perturbations of cTq, over the energy of a a (the superscript 
'o' corresponds to removing one particle from a a and '•' to adding one particle to a a ). Since 
a particle can be only added to do (removed from <7i), in Eq. for / (/1) the first (second) 
term in In is set equal to 1. 

The free energy / of the model can be readily evaluated from Eqs. fl4j) and (jSJ). In 
addition, the transition point fi t can be also obtained — it is the solution of the equation 

III. THE COVERAGE DEPENDENCE OF D c 

Relations (0Q) and (j2J) allows one to calculate an approximate value of the diffusion coeffi- 
cient D c in the local equilibrium limit from the free energy /. Indeed, it suffices to find the 
derivatives 

R = - d -l v = -^ f= (M = ^L m 

dp' X ? N de b K) 



Hence, in combination with Eq. @ that holds analogously also for derivatives of / 18], 
the dependence of D c on \i immediately follows. Consequently, the coverage dependence of 
D c follows upon obtaining the inverse to #(//) and substituting it into D c (/j,), which can be 
easily carried out numerically. 

However, we can also derive explicit finite-size formulas for D c (9), starting from the /i 
dependences of 9, £, and % as yielded by Eqs. (jlj) and ([7]). Without loss of generality, we 



will assume that phase p- (p + ) is stable for p below (above) p t . Then the coverage jump at 
the transition is A9 = 9 + — 6L > 0, where 6± are the single-phase coverages, —df±(p) /dp, 
at fit . 

Notation. As a rule, we use Aq to denote the difference, q + — g_, of single-phase quantities 
q + and q_. 

It turns out that three regimes may be distinguished according to the relative importance 
of phase p_ and p + as given by the weights 

A ±M = l + pile^A/^TV' ( 8 ) 

where the ratio p = V-/v + . Note that < A± < 1 and A_ + A + = 1. In addition, since the 
stable phase has the lowest specific free energy, for p below (above) p t the difference Af 
is positive (negative) and the weight A_ (A+) approaches 1 exponentially fast. At p t the 
specific free energies are identical, and the two weights are equal to 1/2 if p = 1 and to 1/4 
and 3/4 if p = 1/3 or p = 3. 



A. Two-phase regime 

The two-phase regime occurs when 9(p), £(jtf), and x(aO do not reduce to their single- 
phase values, which is true if both A_(/i) and X+(p) are of order larger than A -1 [191 ] . This 
can hold only near the transition — for 

p\ti-m\<6 (9) 

with a small 5 > 0. Indeed, taking the Taylor expansion Af(p) = —A9(p — /it)[l + 0(/3\p — 
p t \)} of Af around p t , h is required that exp(A9NS) is small compared to A. To be specific, 
we shall take 

r 3 InA 

6 = 4A9— (10) 
(instead of 3/4 any number between 1/2 and 1 can be used); then the order of A± is A~ 3//4 
or larger. 

Combining Eqs. (j3J) and (jTJ) with the Taylor expansion 

/V) = /'(^) + ^^(^ - A**) + 0[/3 2 (/i - Ptf] (11) 
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(the prime denotes a derivative with respect to ji or £&), within the two-phase interval (J9]) 

we get 

6(n) = 0-A_(//) + 0+A+(aO + 0(5), 

tfAi) = £_A_(/i)+£ + A + (/i)+0(<5), (12) 
x (/i) = (A0) 2 /3iVA + (//)A_(/.)[l + O(l/ArV4 )]) 

where £± are the single-phase average numbers, df±(n)/deb, of occupied bonds per site at 
/i t . Note that in the two-phase region (jUj) the coverage 6 1 ranges within the interval 

t- < < t+, (13) 

where t± = 9(p t ±5) = 9 ± t A9p ±1 N- 3 / 4 + 0(5). Thus, the coverage attains almost all 
values from to 9 + . 

According to Eq. ffl2|) . the dependence of 9, £, and \ on A 4 i n the two-phase region is 
primarily given by the weights A±. Evaluating A±(0) from the relation 9(p) and equality 
A + + A_ = 1, we readily obtain the coverage dependences of £ and x- Substituting them 
into Eqs. (PJ and (J5J), we arrive at the result 

where P± = 1 — 29 ± + £±/3 is the probability of finding a vacant bond in phase p± at the 
transition and the error term e = 0(N 3 / 4 5) = 0(N~ 1 / 4 \nN). This result holds only for 
coverages in interval ([13]) . 

Formula ( f!4l shows that in the two-phase regime the coverage dependence D c (9) of the 
diffusion coefficient decreases with the system size as 1/N. For a given size, it slowly varies 
if 9 is well between t_ and t + , while it increases as the hyperbole P±/A9N\9 — 0±\ if 9 is 
close to t±. If we neglect thermal effects and the error term e, we can approximately write 



3/26* transition p - P2/3, 

3/2(39 — 1) transition pi/ 3 - pi, (15) 

1 1 9 transition p - p\. 



In Fig. [2] we depict D c (9) in the two-phase interval for four examples of first-order phase 
transitions chosen as follows. 

(Tl) Transition p - p 2 / 3 : > (repulsion), e b = —e t /2 (attraction), and \fi — /i | < £t 
with fi = 3eb/2. 
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FIG. 2: The coverage dependence of the logarithm the diffusion coefficient D c (relative to D®) in 
the two-phase region (|13p for N = 30 x 40. The bond interactions are e b = (—l/2)et < for Tl, 
e b = (-4/3)^ > for T2, e b = (-2/5)e t > for T3, and e h = (-5/3)e t < for T4. The triangle 
interaction f3\et\ = 4 in all cases. Squares correspond to numerical values, whereas dashed lines to 
formula (|14j) . 

(T2) Transition pi/ 3 - p\\ e t < (attraction), Eb = —4^/3 (repulsion), and |/x — /x 1 < \s t \ 
with fi = (9e b + 6e t )/2. 

(T3) Transition Po - Pi'- s t < (attraction), e b = —e t /3 (repulsion), and \fi — /i | < |et| 
with /i = 3e;, + 2e t . 

(T4) Transition p ~ Pi ; £t > (repulsion), = —2^4 (attraction), and |/i — /io| < St with 
/x = 3e b + 2e t . 

The dependence D c (8) is obtained first numerically from Eqs. (PQ), (jl]), and (J7]) and then 
compared to values yielded by formula (|14p with the error term e neglected (in fact, the 
logarithm of D c to base 10 is plotted for better clarity). Obviously, the analytical formula 
very accurately reproduces the numerical results. 

B. Crossover regimes 

At either end the two-phase region is neighbored by a crossover region in which 6(fi), 
£(/i), and xil 1 ) rapidly reduce from two-phase to single-phase values. This corresponds to a 
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decrease in the order of either A_(/i) or A+(/x) from above iV 1 below it. The two crossovers 
take place within the regions 

5 In iV 

S<P\n-to\<d, d= ^§^fi ( 16 ) 

in which the order of A± reduces from TV -3 / 4 to iV~ 5//4 . Using the upper (lower) signs for 
the crossover above (below) p, t , from Eqs. (TjJ and ([7]) we get 

9(ji) = 0± ± x±\^ ~ Ht\ =F Mp ±l e- M ^-^ + e, 

m =e±±C±lM- lh\ T ACp^e-^^ + e, (17) 

X(p) = X± + (A9) 2 j3Np ±1 e~ Ae,3N ^~^ + 0Ne, 

where (± = (d 2 f±/dpd£b) tH , X ± are the single-phase susceptibilities, —d 2 f±(p)/dp 2 , at the 
transition, and the error e = 0(iV~ 3 / 2 ). Thus, in the crossovers the coverage ranges within 
the intervals 

t- < 9 < t-, t + <8 <t + (18) 

with t± = 9{p t i d) — 9± ± (x±/P)d + 0(iV~ 5 / 4 ). Intervals (fl8l) are very narrow and 
concentrated around #±. 

The p dependences of 9, £, and x in Eq. ( TT71) are essentially only via | /x — /x £ | . Evaluating 
the latter [or, more conveniently, evaluating exp(— A9j3\p — fit\N)] as a function of 9, the 
coverage dependences of £ and x can De deduced. Combining them with Eqs. (pP) and (j2J), 
we get 

P± + ££(0 _ ± (£t _ AP)^±M + e 

£> f0) « D°e PfH — — — — — (19) 

with 

u± (9) = W(A9p ±1 C ± e TC±{e - 0±) ) } (20) 

where P'± = —2x± + C±/3 represents the rate of change of P with p in a given phase at the 
transition, the shorthand C± = A9j3N/x±, and W(y) is the Lambert function (the inverse 
to y = WexpW). The upper (lower) signs in formula ( 1191) correspond to the crossover 
around 9 + (0 '_). 

If 9 is near the two-phase region (close to t±), then u± ~ C±\9 — 9±\ ^> 1 so that D c (9) still 
increases as the two-phase hyperbole P±/ A9N\9— 0±\. On the other hand, if 9 is distant from 



the two-phase region (close to r±), then u>± ~ A9p ±1 C± exp(— C±\9 — 9±\) <C 1 so that the 
diffusion coefficient behaves as [P± + {P±/x±){@ ~ ^±)]/(x±//3), i.e., it is linearly disturbed 
from the constant value P±/ (x±/ P)- Thus, within a crossover region the diffusion coefficient 
suddenly changes from the hyperbolic increase to a slight linear increase (or decrease) as 
9 moves from two-phase region across 9± towards a single-phase region. The dependence 
D c (9) in the crossover regions for transitions Tl - T4 is shown in Fig. [3j 



C. Single-phase regimes 

Finally, far from the transition there is a single dominant phase: p- for /i below fit — d 
(or, 9 below r_) and p + for /i above fi t + d (or, 9 above r+). Then #(//), £(//), and x(aO 
reduce to their single-phase values. Indeed, denoting p a the stable phase, Eqs. (TJJ and ([7]) 
yield 

9^) = -^^ + 0(N-^), 

m = d_fM +0{N ^ l (21) 

x(/i) = -^^ + 0(/3Ar-^). 

The coverage dependence of D c in this regime follows by obtaining the inverse to 9(fi), 
using expression (J5]) for / a . This is simple for a = 0, 1, while for a = 1/3,2/3 we may 
take into account that ln(l + x) ~ x because x ~ exp(— const /3) is small. We get the 
approximative formulas 

(1 - 29) /(I - 9f regime of p , 

(5 - 90 + a) (39 - 1 + a)/8g 6 3 a regime of p 1/3 , (22a) 
2/q 3 bV /8q 3 b qf + (2-39y regime of p 2/3 , 

and D c w in the regime of phase pi, where a = [8g^+(l— 36 1 ) 2 ] 1 / 2 and g 6 = exp(— /?£ft). Since 
for phase p\ the approximation of f\ that uses only one-site perturbations is not sufficient 
to get a non- vanishing D c (9), we need to take into account the next dominant contributions 
arising from two-site perturbations (the removal of two particles in a bond). Then an 
additional term (—3/(3) ln[l + exp(— [3 AH*')] appears in fi, where AH'' = 2ji — lleb — I0e t 
is the energy excess of a two-site perturbation of a% over the energy of o\. Applying this 



10 




FIG. 3: The coverage dependence of the logarithm of the diffusion coefficient (relative to D®) in 
the crossover and single-phase regions for the same system size and interparticle interactions as in 
Fig. [2j Squares represent numerical values, whereas dashed and solid lines formula (|14p and (|19p . 
respectively. Dotted lines represent the analytical dependence obtained in the single-phase regions, 
while disks the limiting values. The values A = iV~ 3 / 4 , y± = —1 x 10~ 3 , y% = —7 x 10 -4 , and 
y 4 = 1.2 x HT 11 . 

refined expression for f\, we get 

D° 

D c {0) pa 6 : regime of phase pi (22b) 

qtq^l + 2Aq b q 2 t (l -9) 

with q t = exp(-/3e t ). 

The dependence D c {9) in the single-phase regions for transitions Tl - T4 is detailed 
in the insets in Fig. [31 It turns out that D c {9) does not diverge (as might be incorrectly 
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conjectured from the behavior of D c in the two-phase and crossover regions) but tends to the 
constant values 1, (1 + ^/2q^)/Aq^, 1/ \/2q\q\, and \jq\q\ as 9 approaches 0, 1/3, 2/3, and 1, 
respectively. However, the model does not allow us to analyze the behavior of D c (9) on both 
sides of 9 = 1/3, 2/3 (due to the infinite number of ground states on the lines separating the 
regions of ground states a , a 1 / 3 and cr 2 / 3 , o\). 



IV. CONCLUDING REMARKS 



Figure [3] demonstrates how one analytical approximation quite smoothly takes over from 
another one as the coverage gets from one of the three regions to another. It turns out that 
the agreement increases rather fast with the system size N. 

Although the coverage dependence of D c in the single-phase region as approximated by 
Eq. (1221) is model dependent, the two-phase and crossover regions formulas (tHJ) and (TT9l) are 
actually true for the general case of a first-order phase transition between two phases at low 
temperatures. Indeed, in the local equilibrium limit D c can be approximated by Eq. ([T]) in 
which the factor P is a sum of the probabilities that two or more neighboring lattice sites are 
occupied by adparticles jfj Q, Thus, rather generally, one has P = a + di(df /dpi), 
where a, are constants and Pi are suitable parameters (such as the chemical potential or 
interparticle interactions). Hence, formulas (|14p and f TTTJ]) can be readily extended to this 
general form of P by setting P± = a + Y J i a i(df±/dPi)vt and p '± = Y.i '^ 2 f±/dfidpi)^ t . 
The physical meaning of these quantities remains the same: P± are the single-phase values 
of P and P'± is the rate of change of P with \x in a given phase, both evaluated at the 
transition. 

Finally, note that Eq. (|14p can yield the Langmuir relation, D c (9) = const /(l — 9), if 
phase p + is fully occupied, its factor P + is strictly positive, and the factor P_ in the other 
phase is negligible. This is not the case for the form of P given by Eq. (121). 
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